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The Andreev current through an ultrasmaU NIS junction is calculated in a systematic way by 
means of a nonlinear response approach basing on the elementary Hamiltonian of quasiparticle 
tunneling. The voltage dependence of current and differential conductance as well as the Andreev 
conductance are derived for low- and high-impedance environments, respectively. 
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. Effects of single-charge tunneling at junctions with ultrasmaU capacitances have been studied both theoretically and 
(N ■ experimentally during the last years. For a review see [Q. Modern nanolithography allows the fabrication of tunnel 
: junctions with capacitances less than C < 10~^^F which means that thermal fluctuations can be disregarded at the 
I ' IK scale. Since the tunnel resistances R satisfy the condition R » Rq quantum fluctuations can be neglected. Rq is 
I the quantum resistance Rq = h/e'^. Different types of junctions and several charge transport mechanisms depending 
on whether the electrodes are normal or superconducting ones have been investigated. 

Here the simple example of a voltage driven NIS junction in series with an environment resistance i?^; is studied 
. at zero temperature. At T = there is no charge transport through a NIS junction in terms of quasiparticles for 
\l ' voltages lower than the threshold A/e, where 2 A is the energy gap of the superconductor {A = A{0)). The reason is 
that there are no quasiparticles in the superconductor with energies below A. The quasiparticle current {I)qp reads 



I. INTRODUCTION 



(i),piv) = i^\lv'-(-) e(v--]-[v^^v], (1) 




a • I f A 

g: {i)<ip{v) = - 

^ ■ where the symbol stands for the unit step function. Eq. is valid at T = and in case of a vanishing environment 

Q ' resistance [Re/Rq 0). Since the electromagnetic environment is able to absorb energy it influences the current- 

^ , voltage characteristic of a single junction essentially. In the limit of a high-resistance environment [Re/Rq 00) 

^ the voltage threshold A/e is increased by the Coulomb blockade Ec/e — ej (2C). The corresponding expression of the 

.y—i, , current can be get from Eq. |l] by the replacement V ^ V — Ec/e; 

^: 

^ : {i}„(v) = ji\!{v--) -{-] e{v-=^- = )-[v^-v]. (2) 

But at the moment let us ignore the environment. There is another charge transport mechanism leading to a 
nonvanishing current in the the subgap region (0 < V < A/e). It corresponds to the transfer of two electrons which 
can be converted into a Cooperpair instantaneously. This process, known as Andreev reflection in case of no barrier for 
a long time exists under subgap conditions because it avoids the generation of excited states in the superconductor. 
In this connection the influence of imperfections in the barrier or the electrodes as well as the proximity effect is 
disregarded. The current which is much smaller than ordinary quasiparticle tunneling (Eq. ^ corresponds to a 
process of higher order in perturbation theory and the conductivity emerges to be proportional to R~^. The effect 
was described by means of the Bogoliubov-de Gennes equations as well as in terms of a transition rate formulation 
|^|-^. Here, we are going to show how the nonlinear response approach basing on the elementary Hamiltonian of 
quasiparticle tunneling generates in a systematic and direct way those current contributions (beside others) which 
correspond to the Andreev reflection. After sketching the model (Sec. ^) and basic concepts of nonlinear response 
theory (Se c. P^II[ ), we derive the Andreev current (Sees. IV and ^) for the low-impedance and the high- impedance 



cases (Sec. |VIID, respectively. We restrict ourselves to the concept of ballistic motion of electrons. Finally, a simple 



approximation in case of the high- impedance environment is given (Sec. VII) 
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II. THE MODEL 



The Haniiltonian of a voltage driven NIS junction with environment reads 

H = QV + Hres + Ht , (3) 

where the tunnel Hamiltonian is given by 

Ht = H++H^ , F_ = f| , 

H+ = J2 Tircl^ci^.e'^ . (4) 

ci^a and Cr,CT stand for electron annihilation operators of the left and right electrode, respectively, satisfying anticom- 
mutation relations. The spin is labeled by the subscript a. The sum is taken over momenta I and r and the spin 
(cr =t , i). The reservoir Hamiltonian Hres = Hl + Hn + He consists of terms corresponding to the left {L) and right 
(i?) electrodes and the environment {E) which can be described in standard way • Owing to the phase operator 

tunneling is connected with excitations in the electromagnetic environment. The convention is chosen in such a way 
that a positive voltage favors tunneling from left to right which reduces the junction charge Q by e. The operator H^ 
means tunneling from left to right in contrast to the Hermitian conjugate which describes the reverse process. Tir are 
the tunneling matrix elements. The basic algebra ruling this approach is the relation jl2| 

H± ■ F{Q) = F(Q ±e)-H±, (5) 

where F is an arbitrary function of the junction charge Q. The operator exp(±i<?) changes the macroscopic charge 
Q on the junction by the value =Fe. The algebra corresponds to the elementary commutation relation 

An important point is to take into account the specific feature of a NIS junction. Without loss of generality 
we assume that the right electrode is the superconducting one. Therefore, Bogoliubov transformed operators are 
introduced on the right-hand side in the following way 

4 ,^ = Ur^all^a + Vr^al-r-a , (6) 

where 7^,0- and 7^ ^ are new Fermi operators satisfying common anticommutation relations and the numerical coef- 
ficients Ur^a and are the known BCS coherence factors Q|. We remind of the off-diagonal Hamiltonian of the 
superconducting electrode which becomes diagonal in the new operators 7^ and 7 

Hr = const + ^ Erjl„-fr,a ■ (7) 

r.CT 

The quasiparticlc energies of the superconductor read 



Er = + ^ • (8) 

Hl on the left-hand side is given by 

Hl=Y^ eic\^„ci,a , (9) 

1,17 

where ei are the usual electron energies counted with respect to the Fermi energy. 

III. NONLINEAR RESPONSE THEORY 

The dynamics of our physical system modeled by an unperturbed Hamiltonian H^ and an interaction (tunneling) 
Ht is described by the statistical operator p satisfying the von Neumann equation which corresponds in the interaction 
representation (superscript (J)) to the integral equation 
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p{tY''> -Po-j^ I [H^T\t').p^'\t')W ■ (10) 

It is assumed that the interaction is switched on at t = — oo adiabatically. po is given by the canonical expression 

(/3 = l/{kBT)) 



Po 



The mean value of an arbitrary operator can be calculated in terms of a successive approximation. In case of the 
current operator defined by 



I^^^O^-\[Q,H]^^{H+-H.) (11) 
at in in 



one finds the expansion 



t 

= ! / di'([/«(t),4^)(0])o 



g t t' t" 

i) j dt' j dt' j dt"'([[[/W(i),4'(i')],4^(i")],4^r)])o + ... (12) 

— oo — cxj — oo 

A term of zeroth order is missing because in case of no interaction (tunneling) there is also no current. Now the 
quasiparticle current according to the first order term in Eq. (|l^) can be calculated. The result 

t 

W.P = -^Re j dt'([i/f(i),4^(t')])o (13) 

— oo 

leads at T = to the Eqs. (|l|) and (||), respectively. The Andreev current is contained in the following terms of second 
nonvanishing order 

t t' t" 

(/) = pRe J dt' J dt" J dt"'{{[[[H^^\t),Hi'\t')],HL'\t")],H['\t"')])o 

— OO — OO — oo 

+ m^+\t), H^_'Ht')],H^^\t")],H^_'\t"')])o 

+ m^^\t),H^I\t')],H^J\t")],H^^\t"')])o} , (14) 

where the angle brackets (. . .}o denote averaging with respect to po- In Eq. ( |l^ ) only correlations with vanishing 

signature (+ H and their permutations) are taken into account. The reason is that the separation of the voltage 

dependence in the correlation functions ((. . .)o) with nonvanishing signature (e.g. {H^H^H^H^)o) leads to expres- 
sions containing time dependent (t) terms. Furthermore, a phase remains indeterminated additionally. This is like 



in common Josephson physics where contributions proportional to sin{2eVt/h + <Po) and cos{2eVt/h + <Po) 1 13| arise. 
However, with respect to the Josephson effect the phase is adjusted by current biasing. Therefore, in our opposite 
case of voltage biasing, where the charge becomes adjusted, these terms do not contribute. 

According to the fact that Cooper pairs live in the condensate the corresponding energy balance of the Andreev 
current does not contain any energy contributions belonging to the superconducting bank. It turns out that the 



Andreev current is only determined by the first summand of Eq. (14) because the energy balances of the two other 
summands always depend on the quasiparticle energies Ej.. 
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IV. CALCULATION OF THE CORRELATION FUNCTION 



The voltage dependence can be separated by means of Eq. and in terms of the new time variables 

T = t-t'; t' =t' - t"- t" = t" - t'" 

one gets the Andreev current {lY^^ 

oo oo oo 

(/)(^) = pRe j At j dr' j AT"e~^''^^^+^^'+^"\{T,T'y') (15) 



with 

^(r, r', r") = {[[[H^^'^ (t), H^^\t ~ t)],H^_'^ {t - t - t% H^J\t - r - r' - r")])o ■ (16) 

The hat indicates that the operators only carry the time dependence with respect to the electrodes and the environ- 
ment. The resolution of the interlaced commutators and the separation of the environment lead to 

K = -{H-{t -t-t' - T")H^{t -T- T')H+{t - T)H+{t))o 

X i^e-'^it-r-r'-r")^-^^>{t-r-r')^^^it-r)^^m)^ + 7 further tCrmS . (17) 

Now the operators H± only possess the time dependence with respect to the electrodes 

To simplify matters at first we restrict ourselves to a low- resistance environment {Re ^ Rq)- This means that the 
phase correlation functions can be replaced by 1. The high-resistance case is dealt with in Sec. The calculation of 
the correlation function k is rather lengthy and cannot be shown in detail. The operators H± have to be expressed 
by the elementary operators c, c^, 7 and 7^ which show an exponential time dependence: 

c|t)(t) = e±^-'-*c[5, (18) 
7(t)(t) = e±i^'-*7(t). (19) 

Using Wick's theorem, the correlations can be expressed by the Fermi distribution functions / and the different time 
integrations can be carried out. By means of Dirac's formula, the real part can be taken and the relevant contribution 
coming from the first term in Eq. (|l7| ) turns out to be 

^ J2 fi.ei)f{ev)f{~Er)f{^Er>)urVrUr.Vr,2T*^T;_^Tir.Ti,_r'5{-ei " £/' " 2eV) 
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(20) 



Er -ei~eV Er' - £v - eV Er - ei - eV Er' - ei - eV 
The summation over the spin indices has already been performed. The derivation of Eq. ( po| ) uses the properties 



During the calculation one has to treat lots of terms but only a few of them are nonvanishing and contribute to the 
Andreev current. In similar way as above this analysis can be applied to the other seven terms in Eq. (|l7|). 



V. ANDREEV CURRENT 



Further calculations require an assumption of the momentum dependence of the tunneling matrix elements. With 
respect to quasiparticle tunneling the squared absolute values of tunnel matrix elements iTj^P emerge which are 
approximated by their momentum averages at the Fermi edges of the electrodes (|T/rp) usually. Furthermore, the 
phenomenological tunnel conductance is defined by 
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G 



(21) 



where z/(0) is the density of states at the Fermi level in the normal conducting state. In contrast to this the product 
of tunnel matrix elements in Eq. ( pO| ) is much more complicated because it depends crucially on the nature of the 
electron motion. Some effort has been done to treat this in terms of a diffusive transport ]l^,|7t. However, in case of 
sufficiently small junctions we can assume that the picture of ballistic electron motion is correct which means that 
the scattering of electrons can be neglected. We follow Ref. |^ and approximate the momentum averaged product 
of tunnehng matrix elements 2(T;*T;*_^T;r'T';/_r') in terms of 2{\Tir\'^) ■ 2{\Ti>r'\'^)/N, where iV is a correction factor 
which depends on the geometry of the junction. It corresponds to the number of independent current channels. The 
additional factors 2 are due to the spin. Now, the discrete sums are transformed into integrations where the energy 
densities of states of the superconductor and the normal conductor arise. Finally, the Andreev current is given by the 
following expression: 



(/)(^)(T/) = 



de. 



1 



der J dEr 

— OO 

1 



dErJ{sl)fisl,)fi-Er)f{-E, 
1 1 



Er ~ ei — eV Er 



-A' 



eV Er~ei- eV Er' - £; - eV 



0{\Er\-A)0{\Er,\-A) 
^E^ - Z\2 _ ^2 

S{-ei - ev - 2eV) . (22) 



At zero temperature the Fermi distribution becomes a theta function (/(a;) 
can be carried out. The resulting formula of the Andreev current 



0{~x)) and the energy integrations 



TThG^A 

2 e^N 



log 



1+1 



for 



^ ^ 

— <v < — 

e e 



(23) 



shows logarithmic singularities at ±A/e. This is not surprising because A/e corresponds just to the voltage 

threshold of quasiparticle tunneling and a singularity of such kind is known from inelastic co-tunneling p^ . The 
singularity which is an artefact of perturbation theory shows that still higher order terms are necessary to describe 
the crossover between quasiparticle and Andreev tunneling. The divergence disappears leaving a finite enhancement 
if a finite lifetime broadening is taken into account. This can be done either by hand mm or by performing a 
nonperturbative resummation technique ||l7| , ^ . Of course, the singularity will also be smoothed both due to finite 
temperatures and finite environment resistances. A Taylor expansion in terms of the voltage {V A/e) yields the 
linear behavior 



(/)(^)(F) 



and the Andreev conductance 



^ d(J)(^)(F) 
dV 



nG\ 



nhG^ 



v=o 



2N 



RqG^ 
2N 



(24) 



(25) 



Note, that G*-"^' is proportional to R^^. Fig. ^ shows the Andreev current ( p3| ) and its linear approximation ( |2^ ) 
in the subgap region of quasiparticle tunneling. The scale of quasiparticle tunneling beyond the threshold is about 
several orders of magnitude (factor 2N/{RqG) ^ 1) larger than that of the Andreev tunneling. 



VI. HIGH-IMPEDANCE CASE 



The investigation of the high- impedance case is based on Eq. (17). Now we have to calculate the phase correlation 
functions which were so far replaced by 1. This can be done in Gaussian approximation, for instance by generalizing 
the method presented in Ref. |^ . The phase correlation functions occuring in Eq. (O) depend on the function J which 
contains the information about the structure of the environment. See e. g. Refs. |L0] , pT|Jl9| . In the high-impedance 
case and at T = there is a linear behavior 

J{r) = -z^r (26) 
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which leads to the following expression for the phase correlation function of first term of Eq. ( |l7| ) 



^^_ig.(t-r-r'-r")g-i<i'(t-r-r')gjS>(t-r)gi<J>(t)^^ ^ gi-^ (T+4r'+T") 



Now the first term of Eq. ( p7| ) reads 



8e7r 



J2 f{ei)f{ei')f{-Er)I{-Er')urVrUr>Vr>2TlT;:,_^Tir.Ti.^r'5{-ei - Sy " "leV + 



(27) 



1 



1 



1 



-ei-eV + Ec Er' - ev - eV + E^ Er - ei - eV + E^ Er' - ei - eV + E^ 
Finally, taking all relevant contributions into account, one gets the Andreev current 



(28) 



(/)(-4)(y) = 



dz 



dy 



1 



2{y - w) 



eieV - 2E,) -[V^ -V] , 



(29) 



where the abbreviations w = Ed A and u — eV/ A ~ 2w are used (see Fig. ^, solid line). Eq. ( |29| ) is valid in the 
interval — (Z\ + Ec)/e < V < {A + Ec)/e. Of course, the Coulomb blockade of Andreev tunneling {vj^^^ = 2Ec/e) 
should be smaller than the threshold of quasiparticle tunneling {{A + Ec) / e). This is guaranteed for Ec < A. Though 
the Andreev current cannot be calculated analytically we are able to derive the differential conductance g{V) (see 
Fig. H, sohd line) 



9iV) 



d(/)(^)(V) _ RqG^ 1 



dV 



2N TT^ 



■| + arcsin(?i; — u) ^ + arcsin(ti; + v) 



a/1 - {w + vy 



QieV - 2Ec) + [V ^ -V] 



(30) 



Hence, the Andreev conductance which is defined as the differential conductance at the blockade voltage V^j^'^ turns 
out to be 



RqC^ (l + f arcsinf 



RoG' 



(^ 



arctan ■ 



/l+w 



-9[Vbi )- 2N l-w^ ' 2N l~w-' 



(31) 



The result coincides with that of Ref. Eq. (|2^) is reproduced in the limit Ec — > 0. Fig. (solid line) shows the 

Andreev conductance as a function of w = Ed A. In case that the Coulomb blockade V^i = 2Ec/e tends to the 
voltage threshold of quasiparticle tunneling {A + Ec)/e the Andreev conductance becomes singular [w — > 1). 



VII. APPROXIMATION OF ANDREEV CURRENT 

Since the Andreev current (cf. Eq. (p9|)) is given by an integral which cannot be carried out analytically, we want 
to derive an approximated expression. This can be achieved by taking into account only the relevant singularities of 
the integrand, namely 1/(1 — (w — z)^) and 1/(1 — (t/j + z)^). According to that, the expresion in squared brackets in 
Eq. ( p^ ) is approximated in the following way: 

1 — (w — zY 1 — (w + z)^ ' 

where ai and 02 are constants. Now the integration can be carried out and we get 

TiG'^Aa* fl + w + v l-w + v 
log I • 

e^A^ 47r \l + w — {; 1 — lu — v 

where a new constant a* ~ ai+ 02 is introduced. Furthermore, the factor (1 + w + w)/(l + w — v) in the argument of 
the logarithm can be neglected because it is an irrelevant contribution. The unknown constant a* is determined by 
the demand 



e{v) - [V 



-V] 



(32) 
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approx 



w=0 



hG"^ a* 



w=0 



2N 



which gives a* — 27r^. Hence, the approximated Andreev current reads 

eV 



2N 2e 



A + Ec-eV 



e{eV - 2Ec) -[V ^ -V] 



(33) 



(34) 



We have plotted both the exact and the approximated Andreev current as functions of the voltage for 2Ec < eV < 
A + Ec in Fig. |^. As an example Ec is chosen to be Ec = A/2. The approximated differential conductance is given 
by differentiation 



9iV) 



RqG^ 



1 — w 



— QieV - 2E^) + [V 



'V] 



2N {l-w + v){l-w-v) 

and can be compared with the exact expression (|o|) (see Fig. ||). The approximated Andreev conductance reads 

^approx 



(35) 



(36) 



2N 1-w 

Fig. H shows both the exact and the approximated expressions of the Andreev conductance as a function ofw = Ed A 



ACKNOWLEDGMENTS 



We would like to thank H.-O. MUer for helpful discussions. This work was supported by the Deutsche Forschungs- 
gemeinschaft. 



VIII. CONCLUSION 



In conclusion we have shown how the Andreev current through a NIS junction can be calculated by means of a 
nonlinear response approach basing on the elementary Hamiltonian of quasiparticle tunneling. This method generates 
in a straightforward way all kinds of current contributions including quasiparticle and Andreev tunneling. For low- and 
high-resistance environments we have calculated the voltage dependence of the Andreev current and the differential 
conductance explicitly. Finally, a simple approximation in case of the high-impedance environment has been given. 
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FIG. 1. The Andreev current in the low-impedance case (Eq. ([23[), solid line) and its hnear expansion (Eq. (^, dashed hne) 
are plotted versus the voltage in the interval < eV/A. Current and voltage are counted in units of RqG^ A/{2Ne) and A/e, 
respectively. 

FIG. 2. The Andreev current in the high-impedance case (Eq. (^9|), solid line) and its approximation (Eq. (^), dashed line) 
are plotted versus the voltage in the interval 2Ec < eV < A + Ec for Ec = A/2. Current and voltage are counted in units of 
RQG^A/{2Ne) and A/e, respectively. 

FIG. 3. The differential conductance (Eq. (|30[), solid line) and its approximation (Eq. (|3^), dashed line) are drawn as 
functions of the voltage for 2Ec < eV < A + Ec and Ec = A/2. Conductance and voltage are counted in units of RqG^ /{2N) 
and A/e, respectively. 

FIG. 4. The Andreev conductance (Eq. (pl|), solid line) and its approximation (Eq. (^), dashed line) are plotted as functions 
oi w = Ec/A for < u! < 1. Conductance is counted in units of RqG^ /{2N). 
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